Breakdown of the adiabatic limit in low dimensional gapless systems 
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It is generally believed that a generic system can be reversibly transformed from one state into 
another by sufficiently slow change of parameters. A standard argument favoring this assertion is 
based on a possibility to expand the energy or the entropy of the system into the Taylor series in 
the ramp speed. Here we show that this argumentation is only valid in high enough dimensions 
and can break down in low-dimensional gapless systems. We identify three generic regimes of a 
system response to a slow ramp: (A) mean-field, (B) non-analytic, and (C) non-adiabatic. In the 
last regime the limits of the ramp speed going to zero and the system size going to infinity do not 
commute and the adiabatic process does not exist in the thermodynamic limit. We support our 
results by numerical simulations. Our findings can be relevant to condensed-matter, atomic physics, 
quantum computing, quantum optics, cosmology and others. 



Adiabatic or reversible (also known as quasi- 
stationary) processes by all means play a major role both 
in physics and technology. The adiabatic process is for- 
mally defined as such where no heat is transferred to the 
system from the environment * . Typically such processes 
occur at time scales which are fast enough compared to 
the thermalization times with the environment but yet 
which are sufficiently slow so that the system always 
remains in thermal equilibrium. Another feature of an 
adiabatic process is the conservation of entropy. There 
is a simple general argument showing that slow quasi- 
stationary processes are adiabatic and thus reversible. 
The argument goes as follows: 1 Assume that in an iso- 
lated system some external parameter k is slowly driven 
from some initial value k a to the final one Kg. For sim- 
plicity we assume that k changes linearly in time, though 
it is not necessary, and let 8 be the rate of this change. 
We also refer to S as to the ramp speed. We assume 
that on the way the system does not undergo any dis- 
continuous phase transitions (though second order phase 
transitions are generally allowed) . Then the entropy den- 
sity (or the entropy per unit volume) of the system in the 
final state, Sg, will be a function of this parameter 5 and 
we can expect that for small enough i5 one can expand 
Sb into the Taylor series: 



S B = S A + a'8 + 13' S 2 + 



(1) 



On general grounds we can argue that a' = because the 
entropy can only increase and thus can not be sensitive to 
the sign of 6. Thus at small 6 the excess entropy density 
pumped into the system during this process is 



AS AB « fi'8 2 . 



(2) 



In other words it vanishes quadratically as S — > 0. We 
note there are some subtleties with this expression at zero 
temperature, which we will mention later. 

The adiabatic theorem in thermodynamics is inti- 
mately related to the adiabatic theorem in quantum me- 
chanics, which states that under slow enough external 
perturbations there are no transitions between energy 
levels. Thus, for example, if one starts from a unique 



ground state and adiabatically tunes the system into an- 
other regime of parameters the system remains in the 
ground state and thus no entropy is generated. The 
quantum mechanical adiabatic theorem implies thermo- 
dynamic adiabatic theorem, however, there are some sub- 
tleties involved^. If there is a gap in the system then 
at zero temperature the ground state can be excited 
only through the so called Landau Zener mechanism 3 , 
which gives the exponentially small probability of excita- 
tions and thus the exponentially small entropy increase. 
We note that the fact that this increase is exponential 
rather than quadratic is a peculiarity of the zero tem- 
perature limit. However, systems which have a gap are 
rather exception than the rule. Indeed most of systems 
with broken continuous symmetries have gapless excita- 
tions over the ground state (Goldstone modes). Thus 
solids have phonons (sound waves), ferromagnets and 
anti-ferromagnets have gapless magnon or spin-wave ex- 
citations, superffuids have gapless Bogoliubov excitations 
and so on. Even in gapped systems, for example super- 
conductors, one always has continuum spectrum above 
the gap, which is occupied by quasi-particles at finite 
temperatures. If gapless excitations are present in the 
system then the Landau Zener mechanism does not help 
to protect against creating excitations. However, one can 
generally argue that the available phase space for the 
excitations decreases as the adiabaticity parameter 6 be- 
comes smaller and one can still expect that Eq. @ holds. 

The adiabatic theorem can be also formulated for inte- 
grable systems (the simplest example of these are nonin- 
teracting systems described by harmonic theories) . Then 
the entropy is no longer a good observable because inte- 
grable systems do not thermalize. Instead one can use 
the excitation (quasi-particle) density n cx . The quantum 
adiabatic theorem implies that n ex does not change if 
the process is sufficiently slow. Then on general grounds 
we can expect that the analogue of Eq. {2]| will be 
An ex (S) w (3S 2 . Whether the system is integrable or 
not, its energy is always a good observable, since it is 
defined for any state of the system in or out of equilib- 
rium. In terms of energies the adiabatic theorem can be 
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formulated in the same spirit: 

£b(S) =£ b (0) + (3S 2 , (3) 

where £b(0) is the energy of the state adiabatically con- 
nected to the initial state. 

There is one important caveat in the considerations de- 
scribed above. They are similar in spirit to the mean-field 
argument suggesting existence of the long range order in 
systems with broken continuous symmetries. However, 
we know that at low dimensions either quantum^ or ther- 
mal 5 fluctuations can change the picture and destroy the 
long range order completely. An ultimate reason for this 
is that at low dimensions the density of the low energy 
states p(e) is generically high. Indeed if e(q) cx q z , where 
q is the momentum of an excitation and e(q) is its energy, 
then we necessarily have p(e) cx e d / z ~ 1 , where d is the di- 
mensionality of the system. The quantum or thermal 
fluctuations force excitations to occupy these low energy 
states and if their density is sufficiently high, they can 
qualitatively change the nature of the equilibrium state. 
The same arguments can be applied to the adiabatic pro- 
cess we are interested in. Indeed no matter how slow the 
ramp is (in a gapless system) there will be always created 
some low energy excitations. In low dimensions these ex- 
citations can significantly alter behavior of the system 
and in particular invalidate Eqs. (J2J [3]) . In the remainder 
of the paper we will address this issue in detail. 

The main conclusion of our work is that there are three 
possible regimes of a system response to a slow process: 
(A) mean-field, (B) non-analytic, and (C) non-adiabatic. 
In the first regime (A) we find that indeed the energy of 
the system and other thermodynamic quantities are ana- 
lytic functions of the adiabaticity parameter S and Eq. ([3]) 
is valid. This regime is always realized at sufficiently high 
dimensions. In the second regime (B) one still has the 
adiabatic limit in the sense that there is a well defined 
limit 6 — > but the correction in S to the energy density 
or other quantities behaves nonanalytically: 

£ B (5) ~ £b(0) + P\S\ U , (4) 

where v < 2 is some power which depends on the dimen- 
sionality and on the universal details of the spectrum. 
And finally the third (C) non-adiabatic regime is most 
unusual. There instead of Eq. ^ we have 

£ B (8)^£ B {Q) + (3\5\ U L\ (5) 

where L is the system size and n > is another expo- 
nent. In this regime there is no adiabatic process in the 
thermodynamic limit. In other words the limits <5 — > 
and L — > oo do not commute. This is a striking conclu- 
sion, which states that no matter how slowly we try to 
drive the system, if the latter is sufficiently large, we will 
never reach the adiabatic limit. As we show both analyt- 
ically and numerically such regime is realized in (but not 
limited to) the situation where one starts from an ensem- 
ble of noninteracting bosons in one and two dimensions 



at finite temperature and slowly increases the interac- 
tion strength. We point that existence of the regime C 
is consistent with recent theoretical^ and experimental 8 
works suggesting that there is no adiabatic limit in a par- 
ticular problem of a BCS-BEC crossover in cold atoms. 

Let us make another very important remark. A gen- 
eral argument predicting entropy conservation in a slow 
process assumes that the system is always in thermal 
equilibrium. If, on the other hand, we are dealing with 
an integrable or a nearly integrable system then thermal- 
ization times can be either infinite or very long. However, 
as we argued above the non-integrability of the system 
and its thermalization are not necessary conditions for 
formulation of the adiabatic theorem in the generalized 
sense. The only genuine requirement is that the system 
is initially in a stable equilibrium. In this work we rig- 
orously study the limit where the system either does not 
thermalize or it rethermalizes only after the ramp and 
show that all three regimes are possible. We support our 
results by performing numerical simulations for a partic- 
ular interacting non-integrable system and show that the 
effects of weak nonintegrability do not affect our conclu- 
sions. 

An alternative point of view on the existence of the 
regimes B and C is the breakdown of the linear response 
theory for slow perturbations 9 . In this work we concen- 
trated on the spatially uniform case. However, one can 
anticipate that similar breakdown of the linear response 
will be relevant to a more general class of nonuniform 
low- frequency perturbations. We are leaving the detailed 
analysis of such possibilities for a future work. 

Apart from many applications in condensed matter 
and atomic physics our findings may be relevant to such 
diverse fields as quantum optics, quantum computing, 
in particular for adiabatic quantum computation^; in- 
flationary cosmology, which has the phase of adiabatic 
evolution of Gaussian scale-invariant quantum fluctua- 
tions in the slow-roll approximation^; cosmic microwave 
background radiation^, and Hawking radiation of black 
holes 13 . In all these fields the assumption of adiabaticity 
in slow processes is usually taken for granted, while our 
results suggest that this might not be always the case. 

Analytical treatment. In this paper we consider a spe- 
cific low energy quadratic Hamiltonian: 

where <p q and n g are the coordinate and the conjugate 
momentum. We note that this Hamiltonian describes a 
very wide class of gapless systems. Thus in solids <f> q 
represents the phonon or the plasmon field, in ferromag- 
nets and anti-ferromagnets <f> q describes magnons or spin- 
waves, in superfluids <j> q describes Bogoliubov's excita- 
tions. This list can be easily extended further. Depend- 
ing on the system, the couplings p s and K q have different 
meanings. For example, for superfluids p s denotes the su- 
perfluid density, and K q is related to the compressibility. 
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Dependence of n q on q also varies for different systems. 
Thus n q = const(g) in solids, anti-ferromagnets and su- 
perfluids and n q oc q 2 is ferromagnets and noninteracting 
Bose systems. The reason why the Hamiltonian ((HJ) is so 
generic and applicable to such different situations is that 
it describes Goldstone modes of systems with a broken 
continuous symmetry. We further choose n q = re + Xq 2 . 
This choice allows us to cover all the situations mentioned 
above. In a superfluid re stands for the compressibility of 
the system and we will use this terminology in the paper. 
Let us now imagine that one slowly increases re in time 



re(i) = re + St, 



(7) 



where re is the initial value of the compressibility. In 
principle, one can consider ramps of other parameters, 
but we do not expect any significant changes in the over- 
all picture. In this paper we will focus on positive sign of 
<5, i. e. the situation where compressibility increases in 
time. It can be shown that the opposite process, where 
re decreases, gives similar results up to unimportant nu- 
merical prefactors. We also comment that the choice 
K q (t) = n(t) + \q 2 allows us to analyze interesting pos- 
sibilities from the field theory point of view. Indeed for- 
mally Xq 2 is an irrelevant coupling, which is unimportant 
at low energies (this simply follows from the fact that it 
scales to zero at small q). On the other hand, if one 
starts at zero or very small reo this term dominates the 
behavior of the system at initial times and as we will see 
it qualitatively changes the system's response. Thus A 
plays the role of a dangerously irrelevant variable^ in our 
problem. On the other hand if kq is large then indeed 
one can safely set A to zero. 

Let us first assume that the system is initially prepared 
in the ground state, i.e. consider the zero temperature 
limit. Since the Hamiltonian ^ is noninteracting it does 
not lead to thermalization, therefore as we discussed the 
entropy is not a good concept. One can instead describe 
the degree of non-adiabaticity of the system by comput- 
ing either the density of excitations n cx created during 
the ramp or the energy density (energy per unit volume) 
pumped to the system £ . Since our model is noninteract- 
ing, the evolution of the wave function (or more generally 
density matrix) can be explicitly obtained. We give de- 
tails of such analysis in the section "Methods" . Here we 
will outline the main results. 

Let us first assume that the system is initially prepared 
at zero temperature. Then its wave function factorizes 
to the product of Gaussians corresponding to the ground 
state of a harmonic oscillator for each momentum q (see 
Eq. HH])). We find there are two different regimes of the 
system behavior depending on whether the initial com- 
pressibility reo is finite or zero. The crossover between 
the two regimes occurs at Kq w 8yJ X/ p s . Note that as 



we have 



0. For large initial compressibil- 



ity, Ko S> Kq, the response of the system belongs to the 
A regime in all three spatial dimensions, i.e. the energy 
density behaves as £ oc S 2 . Actually in one dimension 
there is an additional logarithmic correction to this scal- 



ing and A£ oc <5 2 |ln<$|. If reo = then the B regime is 
realized in all three spatial dimensions and the energy 
density scales as 



£ 



(p s A 3 )( d+1 )/ 8 ' 



(8) 



We note that if one analyzes the density of excitations 
n ex then the classification of different regimes is different. 
In particular for large reo one finds that the nonanalytic 
B regime is realized in one spatial dimension and the 
mean-field A regime occurs above two dimensions. For 
reo = the system is in the nonadiabatic C regime for 
d = 1 and it is in B regime in two and three dimensions. 
The fact that the classification of the system's response 
according to the density of excitations is different from 
that based on the energy is not really surprising. Indeed 
primarily low energy modes are occupied and they do not 
contribute much to the total energy. A similar effect also 
occurs in equilibrium. 

If the original system is nonintegrable then it will 
rethermalize at long times. Thermalization in closed non- 
integrable systems is a separate issue, which we are go- 
ing to address in future work. We just point out that 
the Liouville's theorem stating that the entropy of an 
isolated system is conserved in time is generally not an 
issue in non-integrable many-particle systems. Indeed 
the amount of external noise needed to break this en- 
tropy conservation scales to zero exponentially with the 
number of particles so that non-integrable systems are 
never completely isolated. Let us assume that after the 
ramp we can wait long enough so that the system does 
rethermalize and let us see what are the consequences for 
the entropy. Note that since the system is isolated from 
the environment, thermalization occurs at a fixed energy, 
which we just determined above. This energy thus should 
be related to the temperature of the final state Tf via 



£ = C d - Jf 



(K iPs yi 2 ' 



where Cd is a numerical constant. Equating Eqs. 
([9]) we can find that in the case of kq = 



Tf oc ^/KfS 



1/4 A 



3/8 



A 3/8> 



=> S f = A a 



5 d/A 



( P sX 3 ) d / & ' 



(9) 



and 



(10) 



We see that the behavior of the entropy essentially follows 
that of the energy with a slightly different exponent. We 
note that the fact that exponent is different is a peculiar- 
ity of the zero temperature case, where the entropy has 
a singular limit. Apart from that the same classification 
of different regimes A, B, and C applies to the entropy 
of the system, its temperature and other thermodynamic 
quantities. 

Quite similarly one can analyze the response of the 
system to the ramp if it is initially prepared at some 
finite temperature T . Because at finite temperatures the 
excitations are present in the system from the beginning, 
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we will be interested in excess quantities pumped to the 
system during the ramp like An GX , AS, and AS. The 
analysis of finite temperature evolution is straightforward 
and we sketch it in Section "Methods" and Appendix [C] 
The results are again strongly sensitive to the initial 
compressibility. Thus if kq is large then in one dimension 
we get 



AS = A 



I 2 

I Ps^o 



(11) 



In turn, if the system rethermalizes after the ramp, this 
results in the nonanalytic correction to the entropy: 



AS*§A£ 
dS 



A- 



,3/2 ' 



(12) 



In this case obviously the nonanalytic B regime is real- 
ized. In two and three dimensions AS cx S 2 (with loga- 
rithmic corrections at d = 2) and thus we have the mean 
field A regime. 

If kq = then in one and two dimensions the energy 
density diverges with the system size: 



AS = A d 



Ap, 



1Ls 1/3 L 7/3 - d 



1/6 



(13) 



and thus 



AS = A' 



2/3 



d (7-3d)d 
§3(3+1) L 3jd+T) 



(14) 



In particular at d = 1 we have AS cx J 1 / 6 ^ 2 / 3 and at 
d = 2: AS cx (SL) 2 ' 9 . So in one and two dimensions the 
behavior of the system is nonadiabatic and the C regime 
is realized. In three dimensions there is no dependence 
on the system size and we find 



AS = A, 



T 



l/4 A 5 /4 



VS and AS = A' 



Ps 



fp~ s \ 



VS. (15) 



So that the system is in the non-analytic B regime. We 
note again that the density of excitation diverges stronger 
than the energy. In particular at finite T and for kq = 
we find An cx cx (J 1 / 3 ^ 10 / 3 ^, i.e. it diverges with the 
system size in all three spatial dimensions. 

Numerical results: application to interacting bosons. 
While the analysis of the previous sections is formally 
exact, it directly applies only to a special case of an in- 
tegrable harmonic system. Most of the real systems are 
nonintegrable. However, at low energies the harmonic 
approximation is usually very good. We already argued 
that whether the system is allowed to re-thermalize at 
long times or not, the qualitative picture does not change. 
If the thcrmalization occurs then a good measure of the 
heating is the entropy generated during the ramp, if not 
then one should look into the generated density of exci- 
tations. Thermalization in a closed system is certainly a 



very interesting and not very well understood problem, 
which however requires a separate analysis and goes be- 
yond the scope of this paper. So we will focus on the 
analysis of the energy (density) pumped to the system, 
AS, since this observable is not sensitive to the details of 
thermalization. 

To perform numerical simulations we choose the Bose- 
Hubbard model on a square lattice described by the 
Hamiltonian: 



H 



bh 



(ij) 



U(t) 



- 1)' 
(16) 



where aj and aj are the annihilation and creation opera- 
tors of bosons on the j-th site, J represents the tunneling 
matrix element and U is the interactions strength. The 
sum in the first term is taken over the nearest neighbor 
pairs. We take J to be time independent and the interac- 
tion increasing in time according to U(t) = Uo taah((5t) . 
So U(t) first increases linearly in time and then saturates 
at some steady state value. For small enough interactions 
Uq -C J no, where no is the mean number of atoms per 
lattice site, in the quadratic approximation the Hamilto- 
nian (fT6f maps to the Bogoliubov Hamiltonian ([6]) with 
p s w 2 J no, k U , and A ps p s /An^. In order to simu- 
late dynamics of the system we employ the semiclassical 
approach developed by one of us 1 ^. In this approach one 
expands the time evolution of the system in the small 
quantum parameter U/Jno- We note for those more fa- 
miliar with the Keldysh technique 1 ^ that this approach 
treats all classical vertexes exactly and expands the evo- 
lution in quantum vertexes. In the leading order in this 
parameter one obtains the so called truncated Wigner 
approximation (TWA) 16 ' 17 , where the classical fields tplj 

and ipj corresponding to operators aj and aj satisfy the 
time dependent Gross-Pitaevskii equations of motion. In 
the next order the classical fields are subject to a single 
quantum jump during the evolution. We find that while 
TWA approximation is adequate at finite temperatures, 
in the zero temperature limit one has to go beyond and 
add the next correction. This finding agrees with a gen- 
eral statement that the semiclassical approximation can 
break down at long time s 14 ' 18 . We present some details 
of the semiclassical approach in Appendix [D] and here 
only will discuss the results. 

First we look into one dimensional systems, since there, 
according to the theory, we expect strongest effects of 
nonadiabaticity. In order to avoid potential complica- 
tions related to strong quantum effects we choose the 
parameters of the system deep in the superfluid regime 
throughout the entire evolution: no = 20, J = l,Uo — 
0.25 so that the semiclassical parameter Uo/Jno ~ 10~ 2 
(see Ref. [Hj]) is very small. Note that despite Uq = 0.25 
is relatively small the product Uono — 5 is larger than 
J = 1 implying that at long times U(t) w Uq the system 
is in so called quantum rotor regime with healing length 
smaller than the lattice spacing. According to our the- 
oretical expectations the system is mostly excited while 
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the healing length remains large, i. e. when U <C Uq. In 
this regime U{t) — U^tanhSt ~ U^St linearly increases 
with t thus we can directly compare numerical results 
with the analytical predictions of the previous section. 



peratures T ~ 0.01J the behavior of AS is dominated by 
the thermal effects. This trend becomes even more ap- 
parent if we analyze dependence of AS on L (see Fig. ^ . 
In agreement with the analytic results AS strongly grows 
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FIG. 1: Dependence of the energy density pumped into the 
system during the ramp AS on the parameter 8 for different 
initial temperatures and a fixed system size L = 128. The 
other parameters are J = 1, no = 20, and Uq = 0.25. The 
interaction increases in time according to U(t) = UotanhSt. 
The top graph is in the normal scale and the bottom graph 
is in the log-log scale. At finite temperatures the dependence 
of AS on S agrees with the prediction of Eq. (|13[l : AS oc 
(the power 1/3 can be read from the slope of the dependence 
AS (8) in the lower graph). Similarly at zero temperature the 
behavior agrees with Eq. AS oc y8. Note also that at 
small ramps the curves at different temperatures are equidis- 
tant indicating that 8S oc T again in accord with Eq. (|13[) . 

In Fig. Q] we show the dependence of the energy per site 
pumped into the system during the ramp as a function of 
the parameter <5 for different temperatures T at a fixed 
system size L — 128. Note that even at very low tern- 
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FIG. 2: Dependence of AS on 8 for different sizes at fixed 
temperature T = 0.02. As we argued in the previous section 
(see Eq. (|13[0 . there is clearly no thermodynamic limit and the 
heating becomes more severe with the growth of the system 
size. 

with the system size (see Eq. (|T3")) ) and clearly there is no 
thermodynamic adiabatic limit. One can also check that 
the dependence of A£ on L agrees with the prediction of 
Eq. (fT3"|) . Although we do not show the graphs but we 
checked that at zero temperature AS is best fitted by y/5 
dependence and is almost insensitive to the system size, 
again in agreement with the analytic results (jHJ). 

We also performed numerical simulations for the two- 
dimensional Bose-Hubbard model. We work with similar 
parameters as in the one-dimensional case, except that 
at much higher temperature T — 0.2 and slightly smaller 
linear sizes. The reason we had to choose higher T is that 
at a given temperature thermal effects in two dimensions 
dominate the system behavior at much larger linear sizes 
than in one dimension. We find that the results are again 
in a good agreement with predictions of Eq. (p~5|) . In 
particular, we find that AS oc S 1 ^ 3 at small 8 and that 
AS slowly grows with the system size consistent with 
AS ex L 1 / 3 . 

Summary and outlook. In this paper we analyzed 
possible breakdown of the standard adiabatic (or quasi- 
stationary) approximation in low-dimensional gapless 
systems. We explicitly analyzed behavior of harmonic 
systems which are described by quadratic phonon-like ex- 
citations. Despite this particular choice, we point again 
that the Hamiltonian describes a wide class of phe- 
nomena such as phonons in solids, magnons in ferro- and 
antifcrromagnets, Bogoliubov excitations in supcrfluids 
and so force. Using both analytical and numerical meth- 
ods we showed that generically one can have three possi- 
ble heating regimes. The first A regime is mean-field like. 
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L=32x32 
L=64x64 
L=128x128 



FIG. 3: Dependence of A£ on S in two-dimensional system for 
two different sizes. The temperature is fixed at T — 0.2, the 
other parameters are the same as in Fig. [2] The asymptotical 
behavior of AS at small 5 agrees with <5 1//3 dependence pre- 
dicted by Eq. (|f 3[) . Also A£ slowly increases with the system 
size consistent with the analytic prediction AS oc L 1 / 3 . 



There one can apply simple arguments showing that var- 
ious thermodynamic observables are analytic functions 
of the ramp speed 6. This regime is typically realized at 
high dimensions. In the second regime B, which we called 
non-analytic, the energy, entropy and other quantities de- 
pend on 6 in a nonanalytic way: £b{&) ~ £b(0) + (3\5\ v , 
where v < 2. The exponent v depends on the universal 
critical exponents characterizing the gapless phase. And 
finally in low dimensions one can have a very unusual 
regime C, which we called nonadibatic. In this regime 
the limits of 5 — > and the system size L — > oo do not 
commute and for example the energy density behaves as 
£b{5) ~ £s(0) + fl\5\"L n . So no matter how slowly one 
ramps the system, in the thermodynamic limit one can 
never reach the adiabatic regime. 

We did not attempt to classify all possible models and 
model situations where different regimes can be realized. 
This is probably a very difficult task. But as a matter 
of principle we proved that all three regimes are possi- 
ble. We point that there is an interesting connection 
between this work and some earlier works where slow 
dynamics across a quantum critical point was studied. 
See Refs. [19,20.21,22]. In particular, using perturbative 
approach, which we also discuss in Appendix [XJ one of us 
showed^- that the density of excitations generated during 
this process scales as 



(17) 



where z and v are the dynamical and correlation length 
exponent characterizing the phase transition. These re- 
sults were later confirmed by exact methods2ii22 . But 
this is nothing but the B regime of heating, which we 
proposed in this paper. And indeed in Ref. :19] it was 



shown that the scaling (fl7|) is valid below some critical 
dimension, where the exponent of S saturates at 2 and we 
are back to the A regime. We emphasize that our gen- 
eral arguments presented here do not exclude a possibil- 
ity of crossing a continuous second order phase transition 
during the ramp. For example in Refs. [l^lIllHUll] , 
which considered crossings of a critical point in various 
integrable spin chains, the scaling (| 17|) was attributed to 
the Kibble- Zurek mechanism--. However, we stress that 
neither the derivation of this scaling nor the analyzed sys- 
tems do not involve any topological defects. These works 
rather showed that the considered problems belong to 
the nonanalytic regime B of non-adiabaticity. The fact 
that Eq. (|17p is correct is the consequence of the validity 
of the perturbation theory. As we showed in this paper 
the perturbative expression can break down due to the 
bunching effect if the excitations near the critical point 
have bosonic nature. Then the scaling (fTT|) will no longer 
be valid. From results of this paper on general grounds 
one can also expect stronger dependence of An cx on 6 
and possible emergence of the non-adiabatic C regime at 
finite temperatures. 

An obvious outcome of our analysis is that one has to 
be very careful with making statements about adiabatic- 
ity in isolated or nearly isolated gapless low-dimensional 
systems. This can be important in many various situa- 
tions ranging from realizing proposals on adiabatic quan- 
tum computation and preparation of interacting systems 
in a given state via slow ramps to inflationary cosmol- 
ogy and black hole radiation. Perhaps cold atoms are 
the systems where our results can be immediately tested 
in experiments. There one has all the necessary experi- 
mental tools like isolation from the environment together 
with high tunability of parameters of the system and the 
possibility to perform tuning in real time^i. However, 
as coherent dynamic becomes more and more important 
in other fields like nano-physics, quantum optics, cavity 
QED, etc., we expect that our results should be relevant 
to more and more potential applications. 



METHODS. 

Zero temperature. The Hamiltonian is quadratic 
and thus the time evolution can be found exactly. We 
first look into the zero temperature case. For quadratic 
Hamiltonians it is well known that if the initial wave func- 
tion is gaussian, it will remain gaussian for an arbitrary 
time dependence of the parameters of the Hamiltonian. 
Since we assume that the system was originally prepared 
in the ground state the initial wave function is 



mM) = n 



exp 



4cr , 9 



(18) 



where o~o, q = 1/ '(2q)y / 'kq^/ p s . If n changes with time, 
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o q acquires time dependence: 
.da, 



(19) 



This is a Riccati equation, which can be explicitly solved 
through Airy functions. We give details of this solution 
in the Appendix [B] 

It is straightforward to check that the number of exci- 
tations per mode q is related to a q via: 



T cff 



eq 



(20) 



where a q s = l/iJ^cr" 1 ) and cr° q is the equilibrium 
(ground state) value of a in the final state. We point that 
at large values of K the ratio a q s /a q q does not depend 
on k (see details in the Appendix [B]) . The asymptotical 
expressions for n q in the limit of large and small q can be 
easily found from Eq. (|2T))) and the solution of Eq. (TK 
presented in Appendix [Bl In particular, 



(21) 



64q 2 p s nfi 



at q\/p s Ko,q ^ ^ an d 



$1/3 



3^(1/3) q v a (p a K$ iq y/* 



(22) 



in the opposite limit, where T(x) stands for the gamma- 
function. We note that the asymptotics (f2Tj) exactly co- 
incides with the result of the perturbation theory, which 
can be obtained from Eq. (|A1[) . 

The total density of excitations n GX can be obtained by 
integrating n q over q. Let us consider two different cases. 
If kq S> §\f\Jp~ s then n ex is dominated by the integral 
of Eq. (|21|) with the cutoff q ~ 6 and we immediately 
recover the perturbative result, i. e. n cx oc \5\ d for d < 
2 and n cx oc S 2 for d > 2 (in two dimensions ?i ox oc 
5 2 | In (5 1). So we recover that again only type A and type 
B regimes are possible. The situation, where one starts 
in the nonintcracting regime is more complicated. Indeed 
Eq. (j2"2")l indicates that the number of excited modes at 
small q diverges in this case as q~ 4 ^ 3 . In two dimensions 
and above this is an integrable divergence and we again 
recover the perturbative result n cx oc S d / 4 . However in 
one dimension this integral diverges at small q and thus 
should be formally cutoff at q ~ 1/L. In this case one 
finds that n ox oc J^L 1 / 3 . This is precisely the regime 
C, which was missing in the perturbative analysis. We 
already highlighted that the reason why the perturbation 
theory fails is the Bose enhancement of the transition 
probabilities due to their bunching tendency. 

Once we know the complete wave function we can, in 
principle, find arbitrary observables like various correla- 
tion functions, response to external probes etc. However, 



one has to realize that finding wave function was possi- 
ble only due to exact integrability of the harmonic theory. 
So one should look into quantities, which are robust to 
effects of weak nonintegrability. One of such quantities, 
which has a particular importance in thermodynamics is 
the energy density of the system: 



T cff 



(27r) d 



eq 



- 1 



(23) 



where Kf is the final compressibility. The integral above 
converges in all dimensions leading to Eq. ©. 

Finite temperature case. One can generalize the previ- 
ous results to the situation where the system is initially 
prepared in the thermal state characterized by some tem- 
perature T. In this case the excitations are present in the 
system from the beginning and we will be interested in 
analyzing their enhancement during the ramp. Instead 
of the wave function we have to deal with the density ma- 
trix, but essentially the derivation is similar to the zero 
temperature case. We present the details of the analysis 
in the Appendix [C] The result of these calculations is 
surprisingly simple: 



where 



eq 



r q = coth 



q eq 



T=0 



qy/KaTqpl 
2T 



(24) 



(25) 



In other words one needs to take zero temperature 
asymptotics for the width of the wave function and mul- 
tiply them by r q . In the zero temperature limit r q = 1 
and we obviously reproduce the results at T = 0. At high 
temperatures T ^> q^/p s K,Q y q we have 



2T 



q^/plKoTq 



(26) 



and thus a q s diverges at q — > much faster than in the 
zero temperature case. Let us observe that the number 
of excitations of a g-mode is 



T cff 



eq 



r q n q It=0 ^~ 2 ^ Tq 



1). (27) 



The last term is nothing but the equilibrium occupation 
of the g-mode at the initial time. Therefore we see that 
multiplying the zero temperature result for n q by r q we 
immediately obtain the number of additional excitations 
in the g-mode An q generated during the ramp. 

Using Eqs. (JUJ) and 



at qJ PsKq „ > S we get 



A; 



5 2 T 



32q 3 pt 



33/2 7/2 



(28) 



V 0, q 
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and in the opposite limit 



A? 



2n 



<yV3 T 



32/3^(1/3) q ^ p 2 J\ 0q 



(29) 



Now even at finite value of kq the total density of gener- 
ated excitations An GX diverges in one dimension An cx oc 
5 1 / 3 i 1 / 3 and we are in the non-adiabatic C regime. We 
note that formally n cx also diverges at finite tempera- 
ture at d = 1 even in equilibrium, but this divergence 
is much weaker and scales as the logarithm of the sys- 
tem size. In two dimensions we have An cx oc \S\ and 
thus the nonanalytic B regime is realized. And finally 
in three dimensions and above Art ox oc S 2 and thus the 
mean field A regime works (in three dimensions there are 
logarithmic corrections to S 2 scaling). If the initial state 
is noninteracting then the divergence of n q at small q is 
much more severe: n q oc g~ 10 / 3 and in all three spatial 
dimensions the non-adiabatic C regime is realized: 



An cx = A d d — i L i/3+3- f 
p s A 



(30) 



Thus in order to keep the excitation density constant one 
has to scale S as L 3d ~ 10 . This means, for example, that 
in one dimension if the system size increases by a factor 
of two, then in order to keep heating (An ex ) the same 
one has to decrease S by a factor of 128! 
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SUPPLEMENTARY MATERIAL 

APPENDIX A: FERMI GOLDEN RULE 
ANALYSIS 

The easiest way to find the density of excitations n ex 
produced during a slow increase of n(t) in Hamilto- 
nian ^ is to use Fermi golden rule analysis. Indeed, 
since we expect that n cx is small at small S one can ex- 
pect that the perturbation theory in 5 should give a good 
estimate of n ex ((5) at small 5. Using a general expression 
for the density of excitations derived in Ref. [lj| we find: 



1 

32 



d d q 



yexp^-^C 



3/2 



where kq 



(Al) 



Ko + Xq 2 . This expression gives different 
asymptotics in two opposite limits, (i) If 5 p s /A, 



which is the case if one starts in from weak interactions 
kq — > 0, then 



6 d/4 



d/8 x 3d/8 ' 
Ps A 



(A2) 



where A d is a numerical constant. It is easy to check 
that above 8 dimensions the exponent of S saturates at 
2 and does not further change with the dimensionality. 
Expression (|A2I) suggests that in this particular situa- 
tion the nonanalytic regime B is realized in all physical 
dimensions. In one dimension it is particularly hard to 
reach the adiabatic regime since n cx scales only as S 1 ^ 4 . 

We note that the scaling in Eq. (|A2|) is consistent with 
the one obtained in Ref. [19( for the crossing of the sec- 
ond order phase transition: n ox oc S d »/( zl/ + 1 ) j where v 
is the critical exponent characterizing divergence of the 
correlation length. In our case there is a diverging heal- 
ing length £ ~ •*/ X/k instead of the correlation length 
(see Ref. [25[ for details) so that v = 1/2 and given that 
z = 2 in the noninteracting regime one immediately re- 
covers that vf(zv + 1) = 1/4. 

In the opposite limit (ii) where the initial value of k is 
large S <C kq^J p s /X the situation becomes more diverse. 
Thus for dimensions d < 2 Eq. (|A2[) yields 



A' 



d d/2 3(2/2 ' 



(A3) 



On the other hand for d > 2 the exponent saturates and 
we have 



A 1 



-d/2,.d 



Ps 



(A4) 



In two dimensions there is an additional logarithmic cor- 
rection to the scaling (|A3[) . We see that in this situation 
the critical dimension above which the mean field regime 
holds is d* = 2. 

The present analysis can be generalized to other situ- 
ations. For example, in the case of ferromagnets kq = 
and then one can tune A. Then one finds that n cx oc S d ^ 2 
and the critical dimension is d* = 4. We comment that 
one can also consider other scenarios of varying k with 
time. For example, if n oc (5t) r then it is easy to see that 
n GX oc (5 (ir / 2 ( r + 1 ). As r increases the scaling of the density 
of excitations interpolates from (5 d//4 to 5 d ^ 2 and changes 
d* from 8 to 4. 

This simple perturbative analysis shows the existence 
of A (B) regimes for dimensions above (below) some crit- 
ical value d* . This analysis, however, misses the existence 
of the C regime. To justify the validity of the applica- 
tion of the Fermi golden rule one has to require that 
the probability of excitation of each momentum mode is 
small. This requirement breaks down at low energies as 
can be readily seen from Eq. (|A1[) . In the case when 
the excitations have Fermionic character, which is e.g. 
the case for crossing the critical point in the transverse 
field Ising model or the XXZ chain^, the mistake of the 
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perturbative treatment is a simple factor of the order of 
one (see Refs. [T^HOHllIlSl ) - The Goldstone modes de- 
scribed by the Hamiltonian (J6j> on the other hand are 
harmonic oscillators and thus behave as bosons. Bosons 
unlike fcrmions have a bunching tendency, i. e. transi- 
tion probabilities can be significantly enhanced compared 
to the golden rule prediction. 



APPENDIX B: EVOLUTION OF THE WAVE 
FUNCTION UNDER THE RAMP AT ZERO 
TEMPERATURE 



Here we present some details of the solution of Eq. (|19| . 
For convenience we write it here again: 



.da„ 22 1 / % 

= 2 PsQ o- q - -K g (t). 



(Bl) 



This equation can be simplified by first changing inde- 
pendent variable t to n q (t) and then by simple rescaling: 



£1/3 



£1/4 



(B2) 



Under these rescalings we also have k q = k + q 8 ^ 3 . Then 
one can check that Eq. (fT9|) is equivalent to 



.ddq 



(B3) 



This Riccati equation can be explicitly solved in terms of 
Airy functions Ai and Bi: 



. Bi'(-k q ) + a q K\(-k q ) 
Bi(-« g ) + a q Ai(-k q ) 



(B4) 



where a q is an integration constant, which is determined 
from the initial conditions. In the limit k q — > oo ignoring 
unimportant fast oscillating terms we find 



1 



23a„ 



/k q [l 



(B5) 



Note that the real part of l/a q determines \ip\ 2 and thus 
the probability distribution of the corresponding phase. 
The fact that l/cr q — > as k q — > oo should not be surpris- 
ing. Indeed the width of the ground state wave function 
in scaled variables is 



1 



(B6) 



The probability of excitations in the system is determined 
by the ratio of a q and <7 cq , which takes a well defined limit 
at k — ► oo. Introducing a q = l/5ft(<T~ 1 ) we find 



T cir 



cq 



fl« 



23a„ 



(B7) 



The initial condition determining a is: 

r= — _ . Bi(-k , q ) + a q Ai(-ko tq ) 

\f ^0 Q — ^ / \ 7 — ~, T — ■ 

Bi(-Kog) + a q Ai(-k a ^ q ) 
This equation can be inverted to give 

^K ,< z Bi(-K 0i g) - iBi'(-K , 9 ) 



(B8) 



(B9) 



y/k\g Ai(-rtg) - iAi'(-«o, g ) 
In the limit ko^ q <C 1 these equation yields: 

„, R , , 3 2 / 3 r 2 (i/3) _ 

a,«vo + i v K o,g- (BIO) 



Consequently 



T cff 



2tt 



In the opposite limit k 1 ^ 3> 1 one finds a q « i and 

1 



T^ q + 32k 3 



(Bll) 



(B12) 



0. q 



APPENDIX C: EVOLUTION OF THE DENSITY 
MATRIX AT FINITE TEMPERATURES. 



We choose to represent the density matrix correspond- 
ing to the initial thermal state in the Wigner for m 16 ' 28 . 
For the harmonic system described by the Hamilto- 
nian ([6]) one can show that this density matrix factorizes 
into the product of Gaussians: 



W = 



where 



TT— 



■ exp 



4(7, 



o, q r q 



4r„ 



coth 



2T 



, (ci) 



(C2) 



It is well known that in the nonintcracting problem the 
time evolution of the fields <^> q and n q is described by the 
classical equations of motion 2 ^. In particular, 



1 d(f> q 

K q dt 

subject to initial conditions 



d_ 

dt 



+ Ps q 2 = Q 1 



4>q{t = 0) = (f> . q , 4> q (t = 0) = K 0q Tl . 



(C3) 



(C4) 



where 0o,g and Ho, q are randomly distributed according 
to (|Cip . The other important feature of Gaussian en- 
sembles is that in the absence of interactions the Wigner 
distribution (|C1[) always preserves its Gaussian form. 



Therefore finding (</> 2 (i)) and (n 2 (<)) is sufficient to fix 
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the whole distribution function at arbitrary times. Al- 
ternatively one can directly solve the Liouville equation 
for the density matrix in the Wigner form 2 ^ and come to 
the same conclusion. 

A general solution of Eq. (|C3[) is: 



(C5) 



(j>q{kq) = ClAi'(-Kg) + C 2 Bi'(-K g ), 



where as in Appendix [B] we changed the variables from t 
to K q . The integration constants C\ and Ci can be found 
from the initial conditions: 

Ci = 7TK °' 9 ( ?—^- Bi' (-K 0} q) - 7r^o,gBi(-K ,g), (C6) 



v 0,g 



d,K , q 



C 2 = 7r0o,gAi(-Ko, g ) 



g d(p , 

cIkq 



-M'(-k 0iq ). (C7) 



0, q 



From these expressions it is easy to find the asymptotical 
behavior of (0 2 ) at large k and thus find the width of the 



distribution a q s : 



T cff 



2 y/iio,q 



[K ,gBi 2 (-K 0i g) + K , 9 Ai 2 (-K , 9 ) 



+ (Bi'(-« , 9 )) 2 + (Ai'(-^ g )) 2 ].(C8) 

It is straightforward to verify that in the zero temper- 
ature limit Tq — » 1 this expression coincides with the 
previous result (see Eqs. (|B7[) and (|B9[1 V 



FIG. 4: Dependence of the energy density A£ on the S at 
zero temperature with and without the quantum correction 
(|D4[) . For the details of the calculation and the parameters 
of the problem see Fig. [T] Obviously at small values of S the 
TWA breaks down and one has to include the correction (ID4II . 



In the classical limit the bosonic fields ijjj and ipj satisfy 
discrete Gross-Pitaevskii equations: 



Mi 
' dt 



= -jJ2^ + U(t)\^\t 



(Dl) 



ieO-i 



Here the sum in the first term is taken over the nearest 
neighbors of the site j. In the leading order in quan- 
tum fluctuations, so called truncated Wigner approxima- 
tion (TWA), the fields ipj and ip* are subject to ran- 
dom initial conditions, which are distributed according 
to the Wigner transform of the initial density matrix 
W(ip*, ipj). The expectation value of an arbitrary observ- 
able il(a!j, a,j) is given by the average of the correspond- 
ing Weyl symbol (fully symmetrized form of the opera- 
tor) Cl c i(ipt , ijjj) on the solutions of the Gross-Pitaevskii 
equations: 



APPENDIX D: QUANTUM EXPANSION OF THE 
DYNAMICS OF A BOSE-HUBBARD SYSTEM. 



(n(t)) = D^D^W^l^n^it^it)). (D2) 




We find that this level of approximation gives very ac- 
curate results in most of our simulations described in this 
paper. However, at zero temperature case it breaks down 
for very slow ramps and we had to include the next quan- 
tum correction to the TWA. The latter manifests itself in 
the form of a single infinitesimal quantum jump during 
the evolution: 



ipi(t') -> ipi(t') + e x + ie 2 . 



(D3) 



The quantum correction is the evaluated as a nonlinear 
response of f2 c i to such a jumpi^: 
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m))i = - 



o 16 



dei de 2 



def + del 



r 



m 



Numerically both the leading term (Cl(t))o and the 
next correction (f2(t))i are evaluated using Monte-Carlo 
integration schemes. The third order derivatives in 
Eq. (|D4[) are found using finite differences, e. g. 



d 3 n(ei) fi(2ei) - fi(-2ei) - 20(e x ) + 2fi(ei) 
9 3 0(e 1)e2 ) 1 



2,3 



(D5) 



> ,fi(e 1 , e2 ) + 0(e 1 ,-e 2 ) (D6) 



-0(-e 1; e 2 ) - n(-ei, -e 2 ) - 2O( £l ,0) + 2fi(-ei,0) 



It is easy to convince oneself that in order to evaluate 
these finite differences one has to simultaneously solve 
13 Gross-Pitaevskii equations, one for t\ = e 2 = and 
the others for various combinations of ei, e 2 = 0, ±e, ±2e. 
While solving 13 Gross-Pitaevskii equations is certainly 
more time consuming task than solving one, it is still 
tremendously more advantageous than dealing with the 
exact quantum problem. To illustrate the importance of 



quantum correction at zero temperature we show com- 
parison of dependence AS (6) at zero temperature with 
and without this correction (see Fig. 2]). 

Since the initial system is noninteracting, it is straight- 
forward to find the Wigner transform of the density ma- 
trix at finite temperature T. It is more convenient to 
write it in the Fourier space 



-2|-0 9 | 2 tanh 



fo(g) - V 
2T 



(D7) 

where ipk is the discrete Fourier transform of ipj , Z is the 
normalization constant, eo(q) — —JY]j e lqj is the excita- 
tion energy of the Bose-Hubbard Hamiltonian (|16p in the 
absence of interactions and the summation is taken over 
nearest neighbors of site at the origin, n is the chemical 
potential which enforces mean number of particles per 
site uq. We note that in large systems we consider here, 
there is no difference in time evolution between grand 
canonical and canonical ensembles 29 . 
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